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Abstract 

Two comprehensive approaches are considered for constructing projection-based reduced-order computa- 
tional models for linear dynamical systems. The first one reduces the governing equations written in the 
descriptor form, using a Galerkin or Petrov-Galerkin projection onto a reduced-order basis or pair of them, 
respectively. These bases can be constructed by any preferred method such as the Proper Orthogonal Decom- 
position, Balanced Proper Orthogonal Decomposition, or Moment Matching method. The second approach 
transforms first the governing equations into their non-descriptor form before applying the same projection- 
based model order reduction method. For several structural and coupled fluid-structure dynamical systems, 
it is observed that the first reduction approach leads to reduced-order models that exhibit significantly better 
numerical stability and accuracy properties than those resulting from the second one. These observations 
are reinforced by theoretical results that anticipate and confirm the better stability properties obtained in 
general when reducing the descriptor rather than non-descriptor form of the equations governing a linear 
dynamical system. 

Keywords: Balanced POD; Descriptor form; Galerkin projection; LTI dynamical system; Model reduction; 
Moment matching; Non-descriptor form; Petrov-Galerkin projection; POD; Stability 



1. Introduction 

Linear Time-Invariant (LTI) systems are routinely used to model the dynamic behavior of electrical, 
structural, fluid, thermal, and biological systems, to name only a few. Projection-based Model Order Re- 
duction (MOR) methods reduce the dimensionality of LTI systems by projecting them onto carefully chosen 
Reduced Order Bases (ROBs). For many applications, MOR methods have demonstrated the ability to 
preserve the accuracy of large-scale LTI systems [TJ [21 El SI El El [3 [HI - However, perserving their numerical 
stability properties has proved to be a much greater challenge [51 [TU] . A few techniques have been success- 
fully proposed for stabilizing LTI Reduced-Order Models (ROMs) (TTl [TUJ . However, these are not without 
limitations and do not necessarily lead to best-practices in MOR. On the other hand, this paper reports 
on some observations that point to one simple and frequent source of numerical instability in MOR. It also 
presents theoretical results that support these observations and lead to at least two best-practice guidelines 
for performing MOR. 

More specifically, this paper considers two different forms of LTI systems. The first one, known as the 
descriptor form, corresponds to the original form of an LTI system in which a matrix E that is not necessarily 
the identity pre-multiplies the state time-derivative. The second form, known as the non-descriptor form, 
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results from the pre-multiplication of the descriptor form of the LTI system by E -1 , when E is invertible. 
Whereas both descriptor and non-descriptor forms of a given LTI system are mathematically equivalent, 
their projections onto a given ROB generate two ROMs whose numerical properties are observed in this 
paper to differ. Furthermore, it is also proved in this paper that Galerkin-based MOR methods applied to 
the descriptor form of LTI systems deliver ROMs that are less prone to numerical instability than when 
they are applied to the non-descriptor form of these systems. Therefore, two major contributions of this 
paper are the advocacy of Galerkin-based projections over Petrov-Galerkin counterparts when both types of 
projections are feasible for MOR, and their application to the descriptor rather than non-descriptor form of 
LTI systems. 

To this effect, the remainder of this paper is organized as follows. In Section [2] the linear time-invariant 
equations of interest are introduced. In Section [3j the stability of LTI systems is defined for both descriptor 
and non-descriptor forms. In Section [4j the reduction of both descriptor and non-descriptor forms of an LTI 
system is discussed, together with the numerical stability of the resulting ROMs. In Section [5j theoretical 
results associated with the MOR of structural dynamic and coupled fluid-structure systems are presented 
and illustrated with numerical computations. Finally, conclusions are offered in Section [6] 



2. LTI systems 

In this work, first-order LTI systems of the following form are considered 

E^(i)=Ax(t)+Bu(i) 
y(f) = Cx(i)+Du(t), 

where t denotes time, x 6 R™ the vector of state variables, u 6 K p the vector of inputs and y £ K 9 the 
vector of outputs, and E € R nxn , A € R" x ™, B € R nx P, C e R« xn and D e R qx P the time-invariant 
matrices describing the behavior of the physical system of interest. Throughout this paper, E is assumed to 
be non-singular. In Computational Fluid Dynamics (CFD), this matrix can be, for example, the diagonal 
matrix of control volumes arising from the semi-discretization by the finite volume method of the governing 
Euler or Navier-Stokes equations. In finite element methods for structural dynamics (in state space form) 
and heat convection, it is typically the mass matrix. 

Given an LTI system such as ([I]), the following complex- valued transfer function can be defined 

H:seC — > C(sE- A)" 1 B + D e C qxp . (2) 

When E ^ I„, the LTI system ([I]) is said to be in descriptor form. Since E is assumed here to be 
non-singulaij^[ it can be rewritten as 

^(t)=E- 1 Ax(t)+E- 1 Bu(t) 
y(t) = Cx(*)+Du(i), 

which is known as the non-descriptor form of the above LTI system. This form is popular, for example, in 
CFD, as it allows to weight the residual in a given control volume by the inverse of this volume in order 
to emphasize the regions of the computational domain associated with boundary layers where the mesh 
resolution is typically finer. The transformation from the descriptor to the non-descriptor form is also often 
performed to enable the application of a popular MOR method to an LTI system. In this case, a transfer 
function associated with the system ([3| is 

H : s e C i — ► C(s - E _1 A) _1 E _1 B + D e C qxp . (4) 



1 For the more general case where E is singular and therefore there is no alternative to the descriptor form of an LTI system, 
model reduction is studied in |12| and 1131 
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This transfer function is equal to its counterpart ([2|. 

Throughout the remainder of this paper, both systems (JlJ and ([3| are equipped with the initial condition 
x(0) = x € R n . 

3. Stability 

All High-Dimensional computational Models (HDM) of the form ([TJ or ([3| considered in this paper are 
assumed to be stable in the following sense. 

Definition 1. The system Q (respectively j3|) is said to be asympotically stable if lira x(f) = 0, for 

every solution x(f) of E^(t) = Ax(t) (respectively = E _1 Ax(£)). 

til/ 'I '( 

Remark 1. Since systems (UJ) and pi are equivalent, asymptotic stability for one form is equivalent to 
asymptotic stability for the other. 

Next, a Lyapunov function is defined as follows. 

Definition 2. V : R n — > [0, oo) is a Lyapunov function for system (JTJ) (respectively ([3])) if: 



1. V(x) = if and only if x = 

d 



2. For every solution x(i) of E— - (t) = Ax(t) (respectively ~r(t) = E _1 Ax(i)), the following inequality 



holds 

^V(x(t)) < 0. (5) 

Then, the following theorem holds [14]. 

Theorem 1. If the system |l| (respectively ([3])) has a Lyapunov function V, it is asymptotically stable 



Example 1. If V(x) = x (E PE) x, where P is a symmetric positive definite matrix, the inequality (|5 
can also be written in the case of the descriptor form as 

E T PA + A T PE= Q D , (6) 

where Qd is a symmetric positive definite matrix. 

Similarly, if V(x) = x T Px, where P is a symmetric positive definite matrix, the inequality jH]) can also 
be written in the case of the non-descriptor form as 

PE X A + A T E T P= Q ND , (7) 

where Qnd is also a symmetric positive definite matrix. 

4. Projection-based model order reduction 

Whether applied to the descriptor ([I]) or non-descriptor ^ form of an LTI system, a projection-based 
MOR method generates another LTI system of much smaller dimension k « n. In general, such an MOR 
method operates with two ROBs: 

• A trial ROB G R nxk with full-column rank, introduced to approximate the state vector x(i) as 

x(t) w V fc x fc (t). (8) 

The approximate state vector is then uniquely defined by the vector of generalized coordinates x^ G R fe . 
Substituting the above approximation into the original LTI system yields a non-zero residual r(i) G W 1 . 



• A test ROB Wfc with full-column rank as well, introduced to limit the size of the residual r(t) by 



constraining it to satisfy the orthogonality condition W^r(t) = 0. 
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When Wfe 7^ Vfe, the projection-based MOR method is a Petrov-Galerkin approximation method. When 
Wfe = Vfc, it becomes a Galerkin approximation method. 

The next two sections focus on the form of a projected LTI system originally written in descriptor or 
non-descriptor form. Then, the impact of the form of the expression of the LTI system and the type of the 
projection method on the numerical stability of the constructed ROM is discussed in the context of popular 
methods for constructing ROBs. 

4..1. Reduction of an LTI system written in descriptor form 

Once a trial ROB Vfc has been constructed, the state vector is approximated as in pi). Substituting this 
approximation into the original system gives 

EV fe — k - (t) = AV fc x fc (t) + Bu(i) + r(t) 

at (9) 

y(f)=CV t:Xi .(i)+Du(t). 

Next, forcing the residual r(t) to be orthogonal to the test ROB Wfc leads to the following ROM of the LTI 
system 

(W^EV fc ) d ^(t) = (W T k AV k ) Mt)+WlBu(t) 
y{t) = CV kXk (t) + Du(t). 

The counterpart of the above ROM generated by a Galerkin projection is obtained by setting Wfc = Vfc, 
which gives 

(VlEV fc ) ^(t) = (V^AV fe ) xfc(i) + V T k Bu(t) (n) 

y(t) = CV fe x fc (t) + Du(t). 

In general, the trial and test ROBs are chosen to be bi-orthogonal with respect to a well-chosen operator. 
When operating on a descriptor form, Vfc and Wfc are often chosen so that W^EVfc —l k . In this case, the 
ROM of the LTI system can be written as 

^(*)=AfcXfc(*) + Bfcu(t) (12) 
y(t) = CfcX fe (i)+D fe u(i), 

where 

A fc = W^AVfc e R kxk 
B k =WlBeR kxp 

c fc = cVfc e R qxk 



Dfc =D G 



Denoting the approximated high-dimensional state by x(t) = VfcXfc(t), the following equivalent of the 
original high-dimensional system can be reconstructed 

^ (t) = (V fe A) £ (t) + (Vfc W^B) u(t) 
y(f) = Cx(t)+Du(t). 
This system will be particularly useful when studying the stability of the ROM. 
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4-2. Reduction of an LTI system written in non- descriptor form 

For the LTI system expressed in non-descriptor form, the state approximation Q results in 

V fc ^(f) - (E- x AV fc ) Xfc(t) + E- x Bu(t) + r(i) 
y(^CV t x t (t) + Du(t). 

Enforcing the orthogonality of the residual vector r(t) and the test ROB yields 



(15) 



W^E-^Vfe) x fc (t) + (W^E^B) u(t) 
y(t) = CV fe x fc (t)+Du(t). 



(16) 



Here, it is noted that whereas equations (15) are equivalent to their descriptor form counterparts ([9|, the 
reduced LTI system (16) is not equivalent to its descriptor form counterpart (10). In the case of a Galerkin 
projection, this system becomes 



V.V fc - 



dxk 
~dt 



(t) = (V^E^AVfe) xfc(i) + (V^B) u(t) 



y(t) = CV«(f)+Du(i). 



(17) 



If the trial and test ROBs are chosen to be bi-orthogonal with respect to the identity matrix — that is 
W^Vfe = I*., the resulting ROM becomes 



dXfe 

dt 



(t) = A k x k (t)+B k u(t) 



where 



y(t) = C k x k (t)+U k u(t), 



A k = WjE-^Vfe G M fexfc 
B A = W^E _1 B G R fexp 
C fe = CV fc G R qxk 
T> k = D G M 9Xp . 



(18) 



(19) 



As in the descriptor case, the following equivalent of the original high-dimensional system written in 
non-descriptor form can be reconstructed 



dx , 

~dt 



(t) = (VfcW^E-iA) x(t) + (VfeW^E-^) u(t) 



(20) 



y(i) = Cx(t) + Du(t). 



4.3. Stability of the ROM 

In general, a projection-based MOR method with arbitrary ROBs and does not preserve the 
stability of the LTI system to which it is applied. To illustrate this fact, consider the two LTI systems given 
below 



1 
0.5 



dx 



i 

0.6 



-3.5 

-2 



x(i), 



i dx 

and -(*) 



1 -3.5 
0.6 -2 



x(t). 



(21) 



The first one is written in descriptor form, and the second one in non-descriptor form. Both systems are 

asymptotically stable, but the ROM — - — Xl resulting from the Galerkin projection of either of them using 

dt 

is not stable. 



the ROBs Vi = Wi 
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However, the following theoretical results show that when the HDM operators satisfy certain properties, 
stability can be preserved after projection. 

Theorem 2. If the operator A + A T is stable — that is, all of its eigenvalues have non-positive real parts 



— any ROM of the form (12 1 resulting from the Galerkin projection of an LTI system written in descriptor 
form is asymptotically stable. 

Proof. Let V(x) = x T Ex, and let x denote the solution of the ROM problem (14) originating from the 
descriptor form of the LTI system of interest with u(t) = and = V^. Then, 

dV(x) dS. T ^ dx 

— ~ = — Ex + x T E — 
at at at 

= x T (A T V fc V^E) x + x T (EV fe V^A) x 

= x T (A T V fc V^EV fc )x,+xnV^EV fc V^A)x. (22) 
If V^EV fc = I fe and A T + A is stable, then 
dV(x 



dt 



(x T A T V fc ) x fc + (xJV^A) x = (V fc x fc ) J (A T + A) (V fc x fc ) < 0. (23) 



dx 

Hence, V is a Lyapunov function for — = (V^V^A) x which, using Theorem 1, concludes the proof. 

Remark 1. For this choice of function V, there is no analog of the above theorem for the case of a 
Petrov-Galerkin projection of an LTI system written in descriptor form, as the sign of 

^ = xX(A T W t V^E + EV fe W^A)V fc x, (24) 

depends on the choice of the ROBs and Wfc. Moreover, a counterexample will be presented in Sec- 
tion |5.1.2| to show that a Petrov-Galerkin projection of an HDM written in descriptor form and satisfying 
the properties formulated in Theorem 2 does not preserve stability. 

For an LTI system expressed in the non-descriptor form, the following analog to Theorem 2 can be 
established however. 

Theorem 3. If the operator E _1 A + A T E _T is stable, then any ROM of the form (18) resulting from 



the Galerkin projection of an LTI system written in non-descriptor form is asymptotically stable. 

Proof. The proof is identical to that given for the descriptor case using V(x) = x T x. 

Theorem 2 and Theorem 3 highlight an advantage of the Galerkin projection method over its Petrov- 
Galerkin alternative. 



4-4- Sample methods for constructing a ROB 

Proper orthogonal decomposition methods 
The proper orthogonal decomposition method. The Proper Orthogonal Decomposition (POD) method based 
on snapshots [15] computes a trial ROB by compressing the information contained in solution snapshots 
of the system of interest. For the case of LTI systems, these snapshots can be computed either in the time 
or in the frequency domain. To simplify the notation, only the case of a single input (p = 1) and a single 
output (q = 1) is considered here. This case is easily generalizable to multiple inputs and outputs. 

In the time domain, the snapshots are typically obtained by computing the dynamic response of an LTI 
system for a given forcing input and collecting samples {x(^)}^ 1 apa into a (snapshot) matrix 

X = X(t ls . • • ,t Nm J) = [xfa) . . . x(i Nsnap J] , (25) 

where -/V snaps denotes the number of snapshots. 

In the frequency domain, each complex-valued snapshot [3] is computed as the solution of a frequency 
response problem of the form 

(jw,E - A)x(Wi) = B, (26) 



G 



where Ui denotes the sampling circular frequency of interest. Next, these snapshots are stored in a real- valued 
matrix as follows 



X = X(wi J -..,w JW )=[lte(x(a; 1 )) ... Re (x( Wjv .„ p .)) Im(x( Wl )) ... Im (x( Wjv .„ p .))] . (27) 

Whether collected in the time or frequency domain, the snapshots are independent of the form in which the 
LTI system is written because both descriptor and non-descriptor forms of such a system are mathematically 
equivalent. However, as it will be shown below, the trial ROBs Vfc constructed using these snapshots differ. 

In the case of an LTI system written in the descriptor form, when E is symmetric positive definite, the 
POD method proceeds with computing the eigenvalue decomposition 

X T EX = * d A d * dT , (28) 

where X T EX is a small-size matrix. Next, the above decomposition is truncated to the matrices of first k 
eigenvalues and corresponding eigenvectors and ~Vf is constructed as follows 

Vg = X*gAf* (29) 

Alternatively, can be constructed by first computing the Singular Value Decomposition (SVD) of the 
matrix E^X, retaining the first k left singular vectors and defining = E~2 Yj.. Then, a ROM of the 



form (1 1 2h can be built using a Galerkin projection based on V^. 

For an LTI system written in non-descriptor form, the following eigenvalue decomposition is performed 
instead 

X T X = ^d A nd^nd T ( 30 ) 

and VI™ is constructed as in (29]Q A ROM of the form given in (|18[) is then constructed using a Galerkin 
projection based on V£ . 

In summary, because *bf ^ , ^ V£ d and the two ROMs constructed by applying of the POD 
method to both descriptor and non-descriptor forms of the same LTI system using the same snapshots are 
also different. 

The balanced POD method. Information about the outputs of interest can be included in the construction 
of a ROM. For example, the Balanced Truncation (BT) pQ method constructs trial and test ROBs such that 
the ROM resulting from the Petrov-Galerkin projection is balanced in the sense that the most observable 
and controllable states are retained. This method also preserves the stability of the underlying HDM. 
Unfortunately for very large-scale systems, the computational cost associated with the BT method can be 
prohibitive. For this reason, the Balanced Proper Orthogonal Decomposition (BPOD) method [H [TB] was 
developed. This alternative method constructs at a reasonable CPU cost trial and test ROBs that lead 
to an approximately balanced ROM system, but does not necessarily preserve stability (for example, see 



Section 5.1.2). More specifically, BPOD relies on two sets of snapshots. The first one is known as the 
set of primal snapshots. These are identical to the snapshots collected for the POD method. The second 
set of snapshots is known as the set of dual snapshots. These are constructed using a dual LTI system as 
summarized below. 

For an LTI system written in descriptor form, the dual system is 

E T ^(i) = -A T z(i)-C T v(<) , 
dt w w w (31) 

w(t) = B T z(i) + D T v(t), 



2 The trial ROB V? can be equivalently obtained using the SVD of X 



7 



where z 6 R" is the dual state vector, v € K 9 the vector of inputs, and w S K p the vector of outputs. 
Assume that the snapshots are c olle cted in the frequency domain. In this case, each complex-valued dual 
snapshot z d (u>i) associated with (31) is defined as the solution of the following frequency response problem 

(-ME T - A T )z d (uj t ) = C T . (32) 

As described above for the primal snapshots, the dual snapshots are also stored in a real- valued matrix 
as follows 

Z d = Z<V,-.. ,^ Bnap J= [Re(z d ( Wl )) ... Re(z<W„ ap J) Im(zVi)) ■•■ Im (z <W ap j)] . 

( 33 ) 

Next, the SVD of the small-size matrix Z d EX is performed 

Z dT EX = * d S d * dT , (34) 

and the first k singular values and associated left and right singular vectors and are used to 
construct 

Wg = Z rf ¥g£f 4 , and V£ = X d *£sf 4 . (35) 



Then, a ROM of the form ( 10 1 is obtained using a Petrov-Galcrkin projection based on the ROBs W k and 
V fc . 

Similarly, the dual system of an LTI system expressed in non-descriptor form is given by 

* w -ah-. W -o'v(,) (36) 

w(i) = B T E~ T z(<) + D T v(i), 
and the complex- valued dual snapshots z nd (ui) are defined as the solutions of the problems 

(-MIn - A T E- T )z" d (^) = C T (37) 



for multiple circular frequencies Wj. Note that the comparison of Eq. (37) and Eq. (32 ) reveals that the dual 
snapshots associated with both forms of a given LTI system are not identical, but related as follows 

z"V 2 ) = E T z d (-^), i = 1, • • • , A snaps . (38) 

The dual snapshots are stored in the matrix 

Z nd = z n d((Ji; . . . >WJW ) = [ Re ( z » rf ( Wl )) ... Re (z" d ( WjVsnap J) Im (z" d ( Wl )) ... Im (z« d (w Wmp .))] 

(39) 

and the following SVD is performed 

Truncating the first k singular values and associated left and right singular vectors leads to the definition 
of the matrices ^ k d and <i?JJ d and to the construction of the test and trial ROBs as follows 

_ i _ i 

WJJ d = z nd *ff k d 'E k ld 2 , V k ld = X<fr)! d £' fe ici 2 . (41) 

Again, the Petrov-Galerkin projection based on the above test and trial ROBs leads to a ROM of the 
form ( 16 ). 

Focusing for simplicity on the case where the snapshots are computed in the frequency domain, the 
following theorem shows that, unlike the POD method, the BPOD method constructs the same ROM 
whether applied to the descriptor or non-descriptor form of the LTI system of interest. 

Theorem 4. The ROMs obtained by the application of a Petrov-Galerkin projection based on trial and 
test bases computed by the BPOD method to the descriptor and non-descriptor forms of an LTI system of 
interest are identical if the snapshot frequency sampling {wi}^ aps is the same in both cases. 

Proof. See the Appendix. 
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4-. 4 .2. Moment matching methods 

The j-th moment of a transfer function H(s) at s = sq is defined as 

Vj(so) = (-l) j ^(s ), j = 0, 1, 2, ••• (42) 

For the LTI systems considered in this paper, this moment can be expressed explicitly as 

Vo (s ) = C( SQ E-A)- 1 B + D 

ifc(s„) = j!C( So E - A)-(-'"+ 1 )B, j > 1. ' l ' >) 

It can be shown (for example, see [Tl]) that the first k moments of an LTI system can be matched by a 
projection-based MOR method provided that the trial basis satisfies 

range(Vfc) = span{( So E - A^B, • • • , (s E - A)- fc B} . (44) 

Since the transfer function ^ associated with the descriptor form of an LTI system is equal to its coun- 
terpart Q associated with the non-descriptor form of that system, the ranges of the trial bases associated 



with both forms of the LTI system are identical. The right member of (44) describes a Krylov subspace 
defined by the vector (sqE — A) _1 B and the matrix (soE — A)" 1 . Algorithms for constructing this subspace 
can be found in [21 [17j [18], for example. In particular, the Arnoldi method constructs a trial ROB that 
is orthogonal with respect to the identity. For an LTI system written in the descriptor form, a trial ROB 
satisfying V^EVfc = I& can next be constructed by a post-processing step. Then, a ROM can be built using 
a Galerkin projection based on V^. 



5. Numerical stability: Case studies 



5.1. Structural dynamics 
5.1.1. LTI systems of interest 

Linearized structural dynamics problems are typically modeled by second-order LTI systems of the form 



Kz(t) = Fu(i) 
y(t) = Gz(t) + cf t (t) 



(45) 



•Hu(t), 



where z € K™ denotes the vector of structural displacements, u and y the vectors of input and output 
variables, respectively, M, D and K are the mass, damping, and stiffness matrices, respectively, and F, C, 
G and H are time-invariant matrices describing inputs and outputs. The matrices M, D and K are typically 
Symmetric Positive Definite (SPD). Hence, throughout the remainder of this paper, they are assumed to 
have this property. 



Eq. ( 45 1 above can be recast in first-order form by introducing the state- vector 

r/z 

~dt. 
z 



(46) 



containing both the displacement and velocity degrees of freedom (dof). This leads to the following first-order 
LTI system written in descriptor form 



M 




fix 



dt 

y(t) 



(*) = 



D 
J 



-K 




x(t)- 
[C G]x(f)+Hu(f), 



u(t) 



(47) 
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where J € K nx ™ is an arbitrary non-singular matrix. The above system can also be written in non-descriptor 
form as follows 



fix 
~dt 



(*) = 



I 







x(t) 







u(t) 



y(*)=[C G]x(t) + Hu(t). 



(48) 



In general, J is chosen to be the identity matrix. However, the following result suggests that a much 
better choice can be made in the context of constructing a stable ROM. 

Theorem 5. If in the descriptor form (47 1 of the LTI system of interest J is chosen as J = K, then every 
ROM constructed by the application of a Galerkin projection to the HDM (47) is asymptotically stable. 

Proof. Since the matrix D is SPD, 



(49) 



is stable. Therefore, the direct application of Theorem 2 concludes the proof of this theorem. 

Remark 3. The Lyapunov function introduced in the proof of Theorem 2 corresponds to twice the sum 
of the potential and kinetic energies 



D 


-K" 




D 


K 


T 


2D 





K 





+ 


K 














V(x) = x T Ex 



~ dz~ 


T 


M 







~dz~ 










dt 







K 




~d~i 


z 






z 



dz _ ,dz 

— M — 

dt dt 



z'Kz. 



(50) 



Therefore, the decay in time of this Lyapunov function is related to the decay in time of the total energy of 
the structural system due to damping. 

Remark 4. Theorem 3 cannot be applied here because 



-M _1 D 
I 



rvT^K 




-M _1 D 
I 







-M _1 D — D T M _1 
-KM 1 



is not necessarily stable. For example if M = D = K = [1], the resulting matrix 



-M " 1 K 


-2 -1 
-1 



(51) 



unstable. 



Hence, Theorem 2 and this remark highlight an advantage of the descriptor form of first-order LTI systems 
describing structural dynamics problems. 



5.1.2. Analysis of a mass-damper-spring system 

To illustrate the theoretical results presented above, a simple structural dynamics example is considered 
first. It is of the academic type, and therefore offers the advantage of being easily reproducable by the 
interested reader. 

The structural dynamic system associated with this example problem has a one-dimensional chain struc- 
ture (Figure [l|a)) with 20 blocks. Each of these, except the last one, contains four lumped masses, four 
dampers, and five springs. The last block contains one additional spring where an input u is applied (see 
Figure [l|b)) . Therefore, the corresponding one-dimensional first-order LTI system, referred to here as the 
HDM, has n = 160 dof. The element properties of each block ({nii}f =1 , {ci}j =1 , {fcj}| =1 J, where ^5 = 
except for the block where the input is applied, are given in Table [I] The input u is set to a unit step 
horizontal displacement applied at t = 0, and the response of the system is computed for t £ [0, 100] s. 

Two MOR methods are tested using this example: (a) a Petrov-Galerkin projection method using BPOD 
in the frequency domain, and (b) a Galerkin projection method using Moment Matching by Arnoldi's scheme 
(MMA). Each MOR method is applied to the reduction of both descriptor and non-descriptor forms of the 
HDM and the stability of the resulting ROMs is assessed. 

BPOD is applied here using 30 complex-valued primal and 30 complex-valued dual snapshots. These are 
computed in the frequency domain, at evenly spaced frequencies in the interval uj G [0, 27r] rad/s. Then, 
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Figure 1 Mass-spring-damper system 
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Table 1 Mass-damper-spring system element properties. 
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(a) Stability (b) Accuracy 

Figure 2 Mass-damper-spring system: stability and accuracy of the ROMs constructed with BPOD 



ROMs of dimensions k £ {1, ■ • • , 20} are constructed by Petrov-Galerkin projection. The stability of each 
ROM and the relative error of the solution it produces (relative to that produced by the HDM) are reported 
in Figure[2j The reader can observe that, as can be predicted from Theorem 4, the solutions delivered by the 
ROMs constructed from both descriptor and non-descriptor forms of the HDM are identical. The reader can 
also observe for this example that, as stated in Remark 1, the Petrov-Galerkin projection does not always 
preserve the stability of the HDM, even when the properties of Theorem 2 are satisfied. Figure [3](a) shows 
that for the reduced dimension k — 13, the unstable ROMs naturally lead to large relative errors. On the 
other hand, Figure [3^b) shows that for k — 14, the stable ROMs deliver highly accurate solutions. 

For MMA, the moments are matched at so = 2nj and a Krylov subspace of dimension 20 is generated. 
Then, ROMs of dimension k £ {1, • • • , 20} are constructed by Galerkin projection of both descriptor and non- 
descriptor forms of the HDM. Figure [4] reports on the stability and accuracy exhibited by each constructed 
MMA ROM. As predicted by Theorem 5, all ROMs generated by the Galerkin projection of the descriptor 
form of the HDM are stable. The solutions they deliver are also accurate. On the other hand, most ROMs 
generated by the Galerkin projection of the non-descriptor form of the HDM are found to be unstable. 
Consequently, they deliver inaccurate solutions (see Figure [5]) . 

5.1.3. Analysis of a butterfly MEMS system 

Next, a dynamic finite element model of the butterfly MEMS gyro developed in [19] is considered here. 



This model has 17,361 displacement dof and the equations of interest can be written under the form (45 1. 
Therefore, the corresponding first-order LTI system has n — 34, 722 dof. It is referred to here as the HDM 
of this problem. The matrices defining the governing LTI system are part of the Oberwolfach benchmark 
collection [50]. Step input voltages are applied to the undeformed MEMS at t = and focus is set on 
the response of the dynamical system in the time-interval t <E [0, 3] ms. More specifically, the outputs of 
interest are in this case the three axial displacements of the center of a detection electrode. 400 snapshots 
are computed in the time domain and compressed by the POD method, and a Galerkin projection is applied 
to reduce both descriptor and non-descriptor forms of the HDM and construct a suite of ROMs of dimension 
k £ {1, • • • , 20}. The stability of these ROMs and the relative error of the solutions they produce are reported 
in Figure [6j 

The reader can observe that, as predicted by Theorem 5, all ROMs originating from the reduction of 
the descriptor form of the HDM are stable. They are also very accurate. On the other hand, all ROMs 
originating from the non-descriptor form of the HDM are unstable and therefore inaccurate. In particular, 
the time-histories of the outputs predicted by both ROMs of dimension k = 12 and the HDM of dimension 
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Figure 3 Mass-damper-spring system: time-history of the horizontal displacement of the right-most mass 
predicted using the ROMs constructed with BPOD 
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Figure 4 Mass-damper-spring system: stability and accuracy of the ROMs constructed with MMA 




Figure 5 Mass-damper-spring system: time-history of the horizontal displacement of the right-most mass 
predicted using the ROMs constructed with MMA 
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Figure 6 Butterfly MEMS gyro: stability and accuracy of the ROMs constructed with POD 



15 




x 10" 5 



", HDM 


fl 
1 l 


- 1 

i i ■- ROM (descriptor) 


I, 1 1 
1 , II 
' 1 11 


I ROM (non-descriptor) 


' 1 11 1 

I | r i ( 
i , i i , 


pi 


' . 1 ' , 
i , i i , 

1 i 1 1 i 
1,11, 
A' A i A 'a 'a l 
/V / Xi/x'/x F \ / 


UP;: : .; 


pf | 1 1 i 
i , i i , 

' i ' 1 i 
i ,i i . 

' i 1 1 i 
ii ,i i , 
ii i , 
if i , 
* i 1 ' i " 




\' 1 i 
i , 

>. 



Time(s) x1Q -3 
(c) z-axis displacement 

Figure 7 Butterfly MEMS gyro: time-histories of the x, y, and z displacements of the center of a detection 
electrode predicted using the ROMs of dimension k = 12 constructed with POD 



n = 34, 722 are reported in Figure [7J Whereas the predictions delivered by the ROM originating from the 
descriptor form of the HDM reproduces with great accuracy their counterparts delivered by the HDM, those 
obtained using the ROM originating from the non-descriptor form are inaccurate and unstable. 

5.2. Dynamic fluid- structure interaction 

5.2.1. LTI subsystems and coupled LTI system of interest 

Linearized CFD-based computational methods for the solution of coupled fluid-structure interaction 
problems are rapidly becoming the preferred methods for the nutter analysis of modern aircraft in transonic 
and other nonlinear regimes. Here, the linearized, CFD-based computational framework for fluid-structure 
interaction developed first in |21) is adopted for this purpose. In this framework, a structural subsystem is 



modeled as an LTI subsystem of the form (451; its dimension is denoted by n s . For aeroelastic analysis, 



reduced version of this LTI subsystem is typically constructed using a ROB constituted of the first k s natural 



1G 



mode shapes of the structure and a Galerkin projection. Hence, such a reduced system can be written as 



dt 2 



-(*) 



(52) 



where z ks £ ^ ks is the vector of generalized (modal) coordinates, flk B is the diagonal matrix of natural 
circular frequencies of the structure, Dfe s £ ^k B xk s j s £h e generalized damping matrix, p^ is the free-stream 
pressure, rif the dimension of the fluid LTI subsystem described below, x £ R n/ its state, and Ffc s £ K fe = x ™/ 
the generalized Jacobian of the aerodynamic forces acting on the structure with respect to the fluid state 
vector x. 

For simplicity, the flow is assumed here to be inviscid. Hence, it is modeled by the linearization of the 
Euler equations around a steady-state equilibrium characterized by the free-stream pressure p^ and density 
Poo, among others. This modelization leads to a fluid LTI subsystem which can be written as 



^dx / * 



-Ax(t) 




G ks z ks (t) 



(53) 



where E £ J$™/ Xn / is the diagonal matrix of control volumes of the CFD mesh, A <E M. n f xn f is the Jacobian of 
the fluxes with respect to the fluid state evaluated at the steady-state equilibrium solution, and Gfc s £ W l f xks 
and C ke £ M. n f xks are two fluid-structure coupling matrices modeling the interaction of the fluid subsystem 
with the dynamics of the structure. 

Fluid snapshots are generated in the frequency domain by exciting, for a given set of frequencies, the wall 
boundary of the fluid subsystem by one or several structural modal displacement fields. These snapshots are 
compressed using the POD method to construct a ROB Vfc, £ 
using this ROB as 

x(t) w V fc/ x fc/ (t), 

where Xfc, £ 



cxk, 



, and the fluid state is approximated 

(54) 



When the fluid LTI subsystem is expressed in descriptor form as in (53 1, its Galerkin projection onto 
Vfc, leads to the following ROM of dimension kf 



dx fc/ 
dt 



(*) 




(V£ AV fc/ ) x fe/ (i) + .[^»VZ G k z ka (t) + V£ C» 



dz k 



dt 



-(*)• 



(55) 



When Eq. (53) is recast in non-descriptor form and subsequently projected onto the trial ROB V&,, the 
following alternative ROM is obtained 

dx fc . 



dt 



■(*) 



/Poo (^ r T 
poo v 



E^AVfc, x fe ,(t) 



**.(*) 



vr.E- 1 ^ 



dzi 



dt 



-(t). (56) 



In both cases, the reduced fluid subsystem can be written in the generic form 
dx fe 



dt 



it) 




where Afc, £ 



l k f xk f, G kftk3 



I Poo 
Poo 
akfXk 



A kf x kf (t) 




and C fc/i fc s £ 



Poo 
Poo 
kfXk 



Gk f ,k s Zk 3 (t) 



(57) 



Substituting Eq. ( 54 ) into the reduced structural subsystem ( 52 ) leads to the following structural ROM 

dz ka 



where F& 



dt 2 

• k a xkf 



(t) + T) k 



dt 



-(f) +n% s Z ka (t) = y/Po^Fk„k f ^kf(t), 



(58) 



Re- writing Eq. ( 58 1 in first-order form and combining the result and Eq. ( 57 1 into a single LTI system 



leads to the following coupled fluid-structure linear(ized) ROM of size kt + 2k 

d [dz ks ~ 

dt ~dt 

Zfe, 



_:A fe/ _ 

/PooFfe sj fe / 





c k f ,k s 



Ifc. 



-G kf , ks 

-nl a 
o 



dz ks 
dt 

Zfc s . 



(59) 
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(a) CFD surface grid 



(b) Finite element structural model 



Figure 8 High-dimensional aeroelastic model of a wing-store configuration 



Using the above ROM, the onset of flutter can be established by fixing the free-stream density and 
increasing the free-stream pressure p^ until this ROM becomes unstable. At that point, the free-stream 
pressure reaches a value known as a critical free-stream pressure p^ . This fast approach to flutter analysis 



requires however that the coupled fluid-structure ROM (59) be stable outside the flutter point. In [TU], it 
was shown that this in turn requires that the reduced fluid matrix Afc, be stable. Hence, this application 
highlights the importance of requiring that a projection-based MOR method preserves the stability of the 
LTI system or subsystem to which it is applied. 

5.2.2. Flutter analysis of a wing- store- fuel configuration 

Consider first the aeroelastic wing-store-fuel configuration described in [22 and graphically depicted in 
Figure [8| For a fixed altitude, a flight condition for this configuration is defined here by a pair of free-stream 
Mach number and fuel fill level in the store (or tank). The hydroelastic effects due to the presence of 
fuel inside the tank modify the structural properties of the system and affect its aeroelastic characteristics. 
The fluid and structural HDMs developed in [35] for this aeroelastic configuration have the dimensions 
rif = 689,485 and n s = 6,834, respectively. 

For every flight condition of interest, 44 real- valued fluid snapshots are generated by exciting the wall 
boundary of the structural subsystem by each of its first k s = 4 structural mode shapes at each of six 
equispaced reduced frequencies in the interval [0, 0.0125]. Then, these snapshots are compressed by the POD 
method to construct a suite of fluid ROBs of dimension kf e {1, • • ■ , 40}, and a corresponding suite of fluid 
ROMs of same dimension kf is constructed by the Galerkin projection of both descriptor and non-descriptor 
forms of the fluid LTI subsystem onto these ROBs. 



In all cases, the structural ROM is constructed as in (58) with k s — 4 and re- written in first-order form. 

The first considered flight condition is defined by M x = 0.95 and an empty tank. In this case, Figure[9]ja) 
reports on the stability of the constructed fluid ROMs, and Figure |9jb) on the accuracy they deliver for the 
prediction of the critical pressure. Again, all fluid ROMs originating from the descriptor form of the fluid 
LTI subsystem are found to be stable. On the other hand, the fluid ROMs of dimension kf € {29, 36, 37, 38} 
originating from the non-descriptor form of the fluid LTI subsystem are found to be unstable. Consequently, 
each unstable fluid ROM leads t o an erroneous prediction of the critical pressure using the coupled fluid- 
structure (or aeroelastic) ROM (|59l (see Figure Bib)). In contrast, all aeroelastic ROMs of dimension 



kf > 15 originating from the descriptor form of the fluid LTI subsystem deliver accurate predictions of the 
critical pressure. Similar results were also reported in |23j where a preliminary study of this problem was 
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Figure 9 Wing-store configuration (Moo = 0.95 and empty tank): stability of the fluid ROM as a function of 
its size, and accuracy of the critical free-stream pressure predicted using the corresponding aeroelastic ROM 



first performed. 

The second and third considered flight conditions are defined by — 1.1 and 31% fuel fill level in 



the tank, and = 0.75 and 69% fuel fill level, respectively. Figures 10 and 11 report on the stability 
of the constructed fluid ROMs and accuracy of the corresponding aeroelastic ROMs for these two cases, 
respectively. These figures confirm the trends observed for the first flight condition and lead to similar 
conclusions. 

5.2.3. Aeroelastic analysis of an F/A-18 aircraft configuration 



Next, an aeroelastic HDM of a full F/A-18 configuration with tip missiles is considered (see Figure 12 1. 
Here, the dimension of the fluid HDM is n/ = 3,583,095, and that of the structural HDM is n s — 11,290. 
The latter is reduced by the Galerkin projection on a modal basis with k s — 10 structural mode shapes. 

410 real- valued fluid snapshots are computed in the frequency domain by exciting the wall boundary of 
the aircraft configuration with all k s structural modal displacements at 21 equispaced reduced frequencies 
in the interval [0,0.04]. These snapshots are compressed by the POD method to construct a suite of ROBs 
and two corresponding suites of fluid ROMs of dimension kf G {1, • • ■ ,400}: one by Galerkin projection 
of the descriptor form of the coupled fluid-structure LTI system, and one by Galerkin projection of its 
non-descriptor form. The stability of the constructed fluid ROMs is assessed and the obtained results are 
reported in Figure [l3j Once again, all fluid ROMs originating from the descriptor form of the fluid LTI 
subsystem are stable, but more than half those originating from its non-descriptor form are unstable. 



6. Conclusions 

The majority of model order reduction methods do not preserve the stability properties of the large- 
scale linear time-invariant dynamical systems to which they are applied. Stabilization procedures have been 
recently proposed [111 110] to address this issue. However, these procedures have limitations and do not 
necessarily lead to best-practices in model order reduction. This paper contributes to this topic theoretical 
and numerical experiment results that show that the reduction of the descriptor form of large-scale linear 
time-invariant dynamical systems by Galerkin projection preserves their stability under mild assumptions 
that are satisfied, for example, by structural dynamic systems. This suggests that instability in model order 
reduction by Galerkin projection is often caused by the popular reduction of the non-descriptor form of 
the dynamical system of interest, instead of its descriptor form. This suggestion is furthermore supported 
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Figure 10 Wing-store configuration (Moo = 1.1 and 31% fuel fill level): stability of the fluid ROM as a function 
of its size, and accuracy of the critical free-stream pressure predicted using the corresponding aeroelastic ROM 
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Figure 11 Wing-store configuration (Moo = 0.75 and 69% fuel fill level): stability of the fluid ROM as a function 
of its size, and accuracy of the critical free-stream pressure predicted using the corresponding aeroelastic ROM 
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(a) CFD surface grid (b) FE structural model 

Figure 12 Aeroelastic HDM of a full F18/A configuration 
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Figure 13 F18/A configuration at Moo = 0.99: stability of the fluid ROM as a function of its size 
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in this paper by numerous results from different applications obtained using various popular model order 
reduction techniques such as Proper Orthogonal Decomposition, Balanced Proper Orthogonal Decomposi- 
tion, and Moment Matching. Other theoretical and experimental results presented in this paper also lead 
to recommending performing model order reduction on the descriptor rather than non-descriptor form of 
the governing equations, using preferrably a Galerkin projection except when problem specific requirements 
such as those discussed in [23] prevents it. 



7. Appendix: Proof of Theorem 4 

First, the following lemma is stated and proved. 

Lemma 1. In the descriptor case, the dual snapshots satisfy 

Re (z(-cj)) = Re (z(w)) , Im (z(-cj)) = -Im (z(lo)) 

and as a result, the dual snapshot matrix satisfies 



-Wi, 



-lo n . 



snap J = ZVx,--- ,W Wsnap J 



o P 







(60) 



(61) 



Proof of Lemma 1. For ui = 0, the result is straightforward. Let lo ^ 0. Then, Re (z d (w)) and 
Im (z d (w)) are solutions of 



( W 2 E T + A T E~ T A T ) Im (z d ( W )) = -loC t 

Re(z d H) = -E- T A T Im (z d (uj)) 

X ' LO X ' 

Similarly, Re (z d (— u>)) and Im (z d {— lo)) are solutions of 

( W 2 E T + A T E- T A T ) Im (z d (-u;)) = loC t 



(62) 



Re (z d (-uj)) = E" T A T Im (z d (-Lo)) , 



(63) 



which concludes the proof of this lemma. 



Proof of Theorem 4. From the previous lemma and Eq. ( 38 ) , it follows that 











Hence, the product of the dual and primal snapshot matrices is, in the non-descriptor case, 



j rid 



x = 







Z d EX = 







(64) 



(65) 



The last equality is an SVD of the product Z nd X. Eq. kf)b provides another SVD of the same matrix. The 



SVD is not unique; however, assuming that the fc-th and k + 1-th singular values are distinct, the following 
relationship exists between the components of the two SVDs [2"5] 



k 











*^Q, E? d = Eg, *£ d = *^Q, 



(66) 



where Q G IR fcxfe is a block diagonal orthogonal matrix that satisfies QE^Q T = E^. As a result, Q and E d 
satisfy QE^Q T = Ejf * as well. 



22 



Therefore, the trial ROBs for both descriptor and non-descriptor forms are related by 
Similarly, the test ROBs satisfy 



(67) 



W nd = ^ndyndjvnd 2 



E T Z d 











L iV s , 











d~2 



(68) 



Hence, the equivalent HDM (20) reconstructed from the reduction of the non-descriptor form of the 

iPOD can also be written as 

( / 1 = ( Vj2QQ r Wjf EE -1 A) x(t) + (vgQ T QWjf EE _1 b) u(t) 



governing LTI system using BPOD can also be written as 



dt 

y(t) = Cx(t) + Du(t), 



(69) 



(t) = V£W£ A $(t) + V£W£ B u(i) 



y(t)=Cx(i) + Du(t). 



(70) 



The above system is the equivalent HDM ( 14 1 reconstructed from the reduction of the descriptor form of 



the governing LTI system using BPOD, which concludes the proof of Theorem 4. 
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